On the number of real eigenvalues of products of random matrices 
and an application to quantum entanglement 
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The probability that there are k real eigenvalues for an n dimensional real random matrix is 
known. Here we study this for the case of products of independent random matrices. Relating 
the problem of the probability that the product of two real 2 dimensional random matrices has real 
eigenvalues to an issue of optimal quantum entanglement, this is fully analytically solved. It is shown 
that in 7r/4 fraction of such products the eigenvalues are real. Being greater than the corresponding 
known probability (l/\/2) for a single matrix, it is shown numerically that the probability that all 
eigenvalues of a product of random matrices are real tends to unity as the number of matrices in the 
product increases indefinitely. Some other numerical explorations, including the expected number 
of real eigenvalues is also presented, where an exponential approach of the expected number to the 
dimension of the matrix seems to hold. 
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While the study of the spectra of random matrices has 
been extensive and applications have been too numer- 
ous and varied to state briefly [ij, that of products of 
random matrices is relatively fewer, even though it is 
well-motivated 0-0] • For example, products of random 
matrices could describe Jacobian matrices of chaotic sys- 
tems and the rate of exponential increase of the largest 
eigenvalue gives the Lyapunov exponent. A similar con- 
sideration arises for problems involving disordered sys- 
tems, with the matrices being transfer operators and the 
Lyapunov exponent being localization lengths. 

The Ginibre ensemble of random matrices is the sim- 
plest to construct, as these are n x n matrices with all 
the entries being i.i.d. random variables drawn from a 
normal distribution such as A^(0, 1) with zero mean and 
unit variance [Hj. If the entries are complex, the real 
and imaginary parts are independent random variables. 
For the purposes of this work, attention is restricted to 
the real ensemble. It is known that the eigenvalues of 
such matrices can have a significant fraction of eigen- 
values that are themselves real. Explicit expressions for 
Pn,k, the probability that k eigenvalues are real for a ran- 
dom n X n real matrix, have been found. Although these 
are not simple, there are elegant formulae for En, the 
expected number of real eigenvalues as well as the prob- 
ability that there are exactly n real eigenvalues ■ For 
example it is known that: lim„_^oo En/^/n = \f^^T^ and 

It is interesting therefore to study of the number of real 
eigenvalues of products of random matrices. If there are 
K matrices in the product, let the probability that it has 

It is shown below 



k real eigenvalues be denoted by p„ j^, . 

that P2 2 = ""/^j and is therefore larger than the proba- 
bility that a 2 dimensional random matrix has real eigen- 



values, which is p2,2 ~ ^/V^- Numerical results indicate 
that p"^2 monotonically increases to 1 as increases to 
00, thus the probability that there are real eigenvalues 
increases with the number of matrices in the product. 
Numerical results also indicate identical conclusions for 
matrices of dimensions larger than 2, namely that p\n^n 
is a monotonically increasing function of K and seems to 
tend to unity. The distribution of the matrix elements for 
K > \ are not naturally not independent, but that the 
correlations lead to this is a somewhat surprising result. 
Readers not interested in quantum entanglement may go 
directly to the paragraph following Eq. ([5]) . 

(2) 

One direct application of the result for P2.2 to a prob- 
lem in quantum entanglement Q is to find the fraction 
of real "optimal" states of rank-2. A set of pure 

states of two qubits i = 1, . . . , fc} are C-optimal if 

for any {pi, i = 1, . . . , fc, Y.iPi = 0} one has: 



4=1 



k 

E 

4=1 



(10,) (1) 



*e-mail: arul@physics.iitm.ac.in 



C(-) being the concurrence function l^, a measure 
of entanglement between the two qubits. In general the 
R.H.S. is larger than the L.H.S., the concurrence being a 
convex function, and in this sense the set of states leads to 
optimally entangled mixtures if the equality is satisfied. 

Restricting oneself to the set of states that are real in 
the standard basis, it was shown in [l^ that when fc = 2, 
a large fraction (sa 0.285) of pairs of states were in fact C- 
optimal. The sampling of states is such that each of the 
real states is chosen from a uniform distribution on the 
unit 3-sphere , which arises from simply normalization 
of the four real components. Strong evidence was pro- 
vided that the number 0.285 . . ., obtained initially from 

numerical simulations, was in fact (tt — 2)/4. Below it is 

f 2) 

shown that this is in fact P2 2 ~ ^ j whose evaluation then 
confirms the result. For completeness we state that when 
fc = 3 about 5.12% of triples were C-optimal while it was 
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also shown that there was not even one quadrapulet of 
real states that were so. Therefore the set of complex 
states is necessary for there to be C-optimal states in 
general. For k > 2 though, there does not seem to be a 
direct connection to the problem of products of random 
matrices. 

If and 1^2) are an optimal pair satisfying Eq. ([T]) 
we refer to them below as "co-optimal". Such optimal 
pairs satisfy the following conditions [l2j : 

riir22 > 0, and riir22 - rf^ < 0, 

where nj = {(f>i\(Jy ^ay\(l)j). (2) 

Here ay is one of the Pauli matrices. If is a real 
state of two qubits, the concurrence C i\4'i){4'i\) = kiil- 
Let \4>i) = (|00) + |ll))/\/2, be an maximally entangled 
state, so that C (|0i)(^i|) = 1. What characterizes states 
that are co-optimal with this? If 

|,/,2) = a|00) + 6|01) +c|10) {a,b,c,d)eR (3) 

be such a state, then it must satisfy the above conditions 
with rii = —1, r22 = 2(6c — ad) and ri2 — {a + d)/^/2. 
This leads to the following: 

ad~bc> 0, (a + df - 4(arf - be) > 0, (4) 

which begs to be formulated as conditions on the matrix 
M of coefficients defined as 

being equivalent to the requirement that det(Afi) > 
and Ml has only real eigenvalues. Of course if det(Afi) < 
0, the eigenvalues are anyway real. Thus a state is co- 
optimal with a maximally entangled state if its matrix 
of co-efhcients has a positive determinant, yet has real 
eigenvalues. 

Formulated as above, the fraction of states that are 
co-optimal with a maximally entangled state, which we 
can take without loss of generality to be \(j)i) is closely 
allied to the question of the fraction of 2 x 2 real matri- 
ces that have real eigenvalues. The matrix elements can 
be drawn from a normal i.i.d. random process, such as 
N{0, 1). That this gives us the same answer as sampling 
uniformly from the normalization sphere of \cj)2) is evi- 
dent, as the question of reality of eigenvalues of a matrix 
remains independent of overall multiplication by scalars. 
Thus we get the fraction /tt/a of states that are co-optimal 
with a maximally entangled state to be 

f 111 

/V4=P2,2-2=^-2. (6) 

The —1/2 arises as the fraction I/V2 will also include all 
instances when det(Afi) < 0, which are to be subtracted, 
and det(Afi) is equally likely to be positive or negative. 



Thus about 20.7% of real states are co-optimal with a 
maximally entangled one. 

To generalize the above, consider one state as |0i) — 
coseilOO) + sin6'|ll), < 9 < n/4, and the measure fe of 
states that are co-optimal with it, 9 = 7r/4 being what 
was just discussed. That any one state of the pair can be 
chosen such as this follows from Schmidt decomposition. 
The uniform (Haar) distribution on the normalization 
sphere induces an invariant measure, say ^i{9). Then 
the fraction of pairs of states that are co-optimal is given 
by 

(/>= foK0)d9. (7) 
Jo 

The conditions of co-optimality of \(f>i) and a general 
real two-qubit state 102), with rn — — sin20, r22 = 
2{bc — ad) and ri2 = a cos 9 + bsinO, now translate to 
those on the product: 

as det(M2) > and Af2 has real eigenvalues. Once again 
the equivalence of this problem, with (a, b, c, d) chosen 
from the distribution C S{a^ + b^ + + d^ — 1) to one 
where they are chosen from the i.i.d. A^(0, 1) process 
exists. 

Quite independent of the discussion above, but cquiv- 
alently, one may start with a product of two full 
random matrices, say A1A2, and perform a singu- 
lar value decomposition of Ai to get the product 
OiAiOi^^2=OiAiOjA20iOf = OiAii20f . Evidently 
the spectrum of the original product is same as that of 
Ai^2- Here Oi and O2 are orthogonal matrices and 
Ai is a diagonal matrix with positive elements, and 
A2 = O2A2O1. Observe that the if the elements of 
a matrix A are i.i.d. A^(0, 1) distributed, those of the 
products OA, and AO, where O is an arbitrary orthogo- 
nal matrix, are also identically distributed. Therefore it 
follows that A2 has elements that are i.i.d. N{Q, 1) dis- 
tributed. Hence one may well begin with the product in 
Eq. ([8]) without any loss of generality. The diagonal ele- 
ments can be so taken, as any arbitrary singular values 
Ai > A2 may be divided and multiplied by i/Af + Af. 
Identifying cos 9 = Ai/-\/Af + Af , again the resultant ad- 
ditional factor will have no impact on the reality of the 
eigenvalues. 

Let pg be the fraction of matrices M2 that have real 
eigenvalues as {a,b,c,d) are taken from iV(0, 1) and 9 
is fixed. This is realized each time the discriminant 
A2 = (a cos 6* + d sin 6')2 - 4 sin 6* cos 6'(ad -bc)>0. This 
is rewritten as A2 — {a cos 9 — dsm9)^ + 2 sin 20 6c > 0, 
which is a condition on the sum of two statistically in- 
dependent quantities. Using the fact that x = {a cos 9 — 
dsm9) is distributed according to A^(0,1) for all 9 en- 
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ables the following form: 



P9 



e ('4^ 



y ' (27r)3/2 ' ^ > 



where 13 — 1/ sin 20. Note that as — 0, /3 — ?> oo, and 
pe — > 1 . Taking the derivative with respect to /3 converts 
the Heaviside step function into a Dirac delta function. 
Effecting a series of simplifications thereafter, including 
using polar coordinates for y and z results in the following 
remarkably simple equation: 



1 



v^siir^ ( 



d/3 AttJq (sin(/) + /3)3/2- (^0) 

Integrating with respect to /3 and incorporating the 
boundary condition at — gives: 



PS 




1 - — 
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k=0 



k\ 



r(| 



(11) 

The integral in this equation does not seem to acquire 
a simple form except when /3 = 1, which corresponds 
to 6* = 7r/4 and gives = 1/a/2, in agreement with 
the known result, stated previously as P2,2- It follows 
that /e = pe — I is the fraction of states that are co- 
optimal with the state cos 6*100) + sin6'|ll), and one can 
now use Eq. ([7]) to find the fraction of co-optimal pairs. 
The required invariant measure is iJ,{9) = 2 cos 20, this 
follows also as a special case of well-known distributions 
for the singular values of random matrices, for example 
see [1J|. This measure is most well suited to express (/) 



as an infinite series as in Eq. (fTT|) . that may be identified 
with generalized hypergeometric functions. 

Equivalently one may use the integral in this equation 
to express Ji^^^ P0ii{9)d9, the probablity that the product 
of two random 2x2 matrices has real eigenvalues as 

'■'^ sinh^"'^(-y/sin 0) 



(2) _ J_ 



sm( 



(12) 



this follows as the integral can be carried out 
in an elementary way, and also from the evaluation 

^(1-1- sin x)/ sin a; = 2tt. The integral in Eq 
does not appear to be in standard tables, nor fully eval- 
uated by mathematical packages, but as indicated from 
previous work it is in fact simply 77^/2. Therefore it seems 
interesting enough to warrant a more complete evalua- 
tion. The expansion of the inverse hyperbolic functions 
enables the intergal to be written as: 
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FIG. 1: The probability that all eigenvalues of a product of 
K random n dimensional matrices are real, based on 100,000 
realizations. 



Both the generalized hypergeometric functions appear- 
ing here are of the Saalschiitz type, the sum of the top 
rows being 1 less than the sum of the bottom. Theorem 
2.4.4 in 15| can be evoked for such functions, and it is 
remarkable that this is precisely the form of the R.H.S. 
of the identity therein, which results in its evaluation as 
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(14) 



Here use is made of an identity of Gauss for 2F1 at argu- 

(2) 

ments of unity [15], and leads to P2 2 — ^/4, and hence 
finally (/) = (n - 2)/4. 

The generalizations dealt, with numerically below, are 
to products of more matrices as well as to higher dimen- 
sions for a fixed number of products. The behavior of 
^2^' for i^T > 2 is seen in Fig ^ and shows this mono- 
tonically increasing with K. In the same figure is also 
shown the corresponding probablities that all the eigen- 
values are real for such products of 3 and 4 dimensional 
matrices. This increase in the probability that all eigen- 
values are real is also reflected in the expected number 
of real eigenvalues. This is shown in Fig. ([5]) where this 



number: = J2k=o ^PlCk plotted as a function of n 

for fixed values of number of products K in the top panel. 
In the bottom panel the expected number is shown as a 
function of K for 2 and 4 dimensional matrices. The 
exponentially fast approach to the dimension of the ma- 
trix is to be noted. Thus n — E^'^ ~ exp(— 7„iir) seems 
to hold with 7„ decreasing with increasing n (72 and 74 
are the slope values in Fig. ©. Numerical results also 
indicate that the probability that there are k real eigen- 
values for k < n while not necessarily monotonic, does 
eventually vanish with the number of products, leaving 
the dominant case as the one with all eigenvalues real. 



This is illustrated in Fig. ^ where the quantity 
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FIG. 2: The expected number of real eigenvalues of products 
of K random n dimensional matrices, based on 100,000 real- 
izations. The top panel shows variation with the n for fixed 
values of K, while the bottom one displays the exponential 
approach of the expected number of real eigenvalues to the 
dimension of the matrix as K increases. 
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FIG. 3: The probability that k eigenvalues of a product of 
K random 8 dimensional matrices are real, based on 100,000 
realizations. The k = Q case is barely seen in this scale. 
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FIG. 4: The eigenvalues of K products of 10 dimensional 
random matrices, after they have been divided by the corre- 
sponding Frobenius norms. The real and imaginary parts are 
plotted for 1000 realizations of such products. 



of if, namely 1, 2, 5, and 10, and the distortion from an 
approximately circular law is evident with the formation 
of two lobes. The eigenvalues are divided by the norm of 
the products of the matrices so that the values are not 
exponentially increasing and can be compared. 

In conclusion and summary, the number of real eigen- 
values for products of real random matrices has been 
studied. The case of products of two 2 dimensional ran- 
dom matrices was fully analytically solved and it was 
shown that in the fraction of 7r/4 cases, the matrices had 
real eigenvalues. This solved a problem of entanglement, 
where it was shown that the fraction of optimal pairs of 
two qubit states is therefore (tt — 4)/2. Generalizations 
show that with increasing number of products, all the 
eigenvalues tend to be real with probability approaching 
unity. Needless to say, the numerical results pose inter- 
esting challenges, as the resulting matrices have highly 
correlated matrix elements. 



plotted for 8 dimensional matrices for k — 0,2,4,6, and 
fc = 8, corresponding the probability of finding k real 
eigenvalues. 

For a fixed dimensionality as the number of products 
increases more and more eigenvalues "condense" from the 
complex plane onto the real axis. The distribution of 
the eigenvalues hence changes significantly as well. For 
a single random matrix the eigenvalues are roughly dis- 
tributed according to the circular law [l^, Il3| , while the 
real eigenvalues are asymptotically uniformly distributed 
0. In Fig. Q are shown the eigenvalues of products of 
10 dimensional real matrices. This is shown for 4 values 
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